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Abstract 


In this paper, we report a version of the Space-Time Conservation Element and So- 
lution Element (CE/SE) Method in which the 2D and 3D unsteady Euler equations are 
simulated using structured or unstructured quadrilateral and hexahedral meshes, respec- 
tively. In the present method, mesh values of flow variables and their spatial derivatives 
are treated as independent unknowns to be solved for. At each mesh point, the value of 
a flow variable is obtained by imposing a flux conservation condition. On the other hand, 
the spatial derivatives are evaluated using a finite-difference/weighted-average procedure. 
Note that the present extension retains many key advantages of the original CE/SE method 
which uses triangular and tetrahedral meshes, respectively, for its 2D and 3D applications. 
These advantages include efficient parallel computing, ease of implementing non-reflecting 
boundary conditions, high-fidelity resolution of shocks and waves, and a genuinely mul- 
tidimensional formulation without using a dimensional-splitting approach. In particular, 
because Riemann solvers — the cornerstones of the Godunov-type upwind schemes, are not 
needed to capture shocks, the computational logic of the present method is considerably 
simpler. To demonstrate the capability of the present method, numerical results are pre- 
sented for several benchmark problems including oblique shock reflection, supersonic flow 
over a wedge, and a 3D detonation flow. 

1. Introduction 

The Space-Time Conservation Element and Solution Element (CE/SE) Method, orig- 
inally proposed by Chang [1-13], is a new numerical framework for solving conservation 
laws. The CE/SE method is not an incremental improvement of a previously existing 
CFD method, and it differs substantially from other well-established methods. The CE/SE 
method has many nontraditional features, including a unified treatment of space and time, 
the introduction of conservation element (CE) and solution element (SE), and a novel 
shock capturing strategy without using Riemann solvers. Note that conservation elements 
are nonoverlapping space-time subdomains introduced such that (i) the computational do- 
main is the union of these subdomains; and (ii) flux conservation can be enforced over 
each of them and also over the union of any combination of them. On the other hand, 
each solution element is a space-time subdomain over which any physical flux vector is 
approximated using simple smooth functions. In general, a conservation element does not 
coincide with a solution element. 

To date, numerous highly accurate CE/SE steady and unsteady solutions with Mach 
numbers ranging from 0.0028 to 10 have been obtained without using preconditioning 
or other special techniques [1-26]. The flow phenomena modeled include traveling and 
interacting shocks, acoustic waves, shedding vortices, detonation waves, and cavitation. In 
particular, the rather unique capability of the CE/SE method to resolve both strong shocks 
and small disturbances (e.g., acoustic waves) simultaneously has been verified through 
several accurate predictions of experimental data [15-17]. Note that, while numerical 
dissipation is required for shock resolution, it may also result in annihilation of small 
disturbances. Thus a solver that can handle both strong shocks and small disturbances 
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simultaneously must be able to overcome this difficulty . The design principles of the 
CE/SE method have been extensively illustrated in the cited references. In this paper, a 
brief description of the CE/SE method is provided as the background of the present work. 

Perhaps, one of the most important features of the CE/SE method is the adoption 
of an integral form of space-time flux conservation as the cornerstone for the subsequent 
numerical discretization. Note that one derives the conventional finite-volume methods 
based on Reynolds’ transport theorem [27], in which space and time are treated separately. 
As will be shown shortly, this separate treatment of space and time imposes a restriction 
on the space-time geometry of finite volumes and, as a result, classical Riemann problems 
arise natually in the course of flux evaluation across an interface. In contrast, due to 
its unified treatment of space and time, Chang’s flux conservation formulation allows a 
choice of the space- time geometry of CEs that render it unnecessary to solve Riemann 
problems. To clarify this fundamental difference, in this Introduction, we will first review 
the conventional integral form for hyperbolic conservation laws in Sec. 1.1 as a contrast 
to Chang’s integral form which is described in Sec. 1.2. The original CE/SE method is 
reviewed in Sec. 1.3, and the objectives and outline of the present work are presented in 
Sec. 1.4. 

1.1. Conventional Finite Volume Methods 

Consider the differential form of a conservation law, i.e., 

! + v.£ = o (i.i) 

where (i) u is the density of the conserved quantity; (ii) h is the spatial flux vector; and 
(iii) V- is the spatial divergence operator. Note that, in order to distingush a spatial object 
from a space- time object (see below), hereafter the former will be denoted by an underline. 
By using the Reynolds’ transport theorem, one can obtain the conventional integral form 
of Eq. (1.1), i.e., 

f u dv + (f> h • ds = 0 (1.2) 

dt Jv Js(V) 

where (i) V is a fixed spatial domain (i.e., a “control volume”); (ii) dv is a spatial volume 
element; (iii) S(V) is the boundary of V; and (iv) ds = dan with da and n, respectively, 
being the area and the unit outward normal vector of a surface element on S(V). By 
integrating Eq. (1.2) over the time interval (t 5 ,ty), one obtains 
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The discretization of Eq. (1.3) is the focus of conventional finite volume methods [27]. 

1.2. The Space-Time Flux Conservation Formulation 
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Let, (i) E\ denote an N-diinensional Euclidean space in which Xi,X 2 , • ■ • , r/v-i are 
spatial coordinates and x ,v = t\ (ii) V- be the divergence operator in £\y ; and (iii) h 
( h, u ). Then Eq. (1.1) implies V • h = 0. As a result, Gauss’ Divergence Theorem in En 
implies 

(p h ■ ds = 0 (1.4) 

Js(V) 

As depicted in Fig. 1, here (i) S(V ) is the boundary of an arbitrary space-time region V in 

E,v, and (ii) ds = dan with da and n, respectively, being the area and the unit outward 

— * 

normal of a surface element on S(V). Note that: (i) because h • ds is the space-time flux 
of h leaving the region V through the surface element ds , Eq. (1.4) simply states that 
the total space-time flux of h leaving V through S(V’) vanishes; and (ii) all mathematical 
operations can be carried out as though Ejv were an ordinary N-dimensional Euclidean 
space. 

Let N = 2. For this case, (i) X\ = x and x 2 = t ; (ii) h = h x ; (iii) V • h = dh x /dx\ and 
(iv) a “surface element” on S(V) and the “area” of this element reduce to a line segment 
and the length of this segment, respectively (see Fig. 1). Note that, for an arbitrary V, 
the spatial projection V(t) of the cross-section of V at time t generally varies with t. The 
exception occurs only if V is a cylinder with its axis being parallel to the time axis, such 
as the rectangle ABCD depicted in Fig. 2(a). In this case, V_{t) is independent of t and 
thus it can be considered as a “control volume.” 

Let V be the rectangle ABCD depicted in Fig. 2(a). Then S(V) is formed by the line 
segments AB , BC , CD and DA. Let (i) t = t 3 at CD; (ii) t = tf at AB; (iii) x = x 9 at 
BC ; and (iv) x = x/ at DA. Then because h = (h x ,u), with the aid of Fig. 2(a), Eq. (1.4) 
implies 



Note that Eq. (1.3) reduces to Eq. (1.5) for the ID unsteady case in which (i) V is the 
spatial cylinder of constant cross-section depicted in Fig. 2(b); (ii) u = u(x,t)-, and (iii) 
h = (/ij.,0,0) with h x = h x (x,t)- 

Note that generally the discretization of Eq. (1.3) is carried out by dividing the entire 
space-time computational domain into space-time CEs. Each CE is a cylinder in space-time 
with (i) its spatial projection being the control volume V_, and (ii) its top and bottom faces 
representing two constant time levels. Because the control volume is a fixed spatial domain, 
these CEs generally are stacked up exactly on the top of each other, i.e., no staggering of 
CEs in time is allowed (see Fig. 2(c) for the N = 2 case). With this arrangement of CEs, 
the vertical interface that separates any two neighboring columns of CEs will always be 
sandwiched between two neighboring columns of mesh points (marked by dots in Fig. 2(c)). 
As such, flux at the vertical interface of two neighboring CEs generally must be evaluated 
by interpolating the data from these two CEs. How this interpolation should be carried 
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out properly under varying solution behavior is a difficult problem. As will be shown 
immediately, with a new .space-time arrangement of CEs and mesh points, and a proper 
definition of SEs. the above difficult interpolation problem can be bypassed completely. 

1.3. The CE/SE Method 

As an example, the CE/SE method will be described by considering the PDE 


du d(au ) 

— + — = 0 

dt dx 


( 1 . 6 ) 


where a is a constant. Obviously the integral form of Eq. (1.6) is Eq. (1.4) with N = 2 
and h. = (au, u). 

To proceed, let 'I' denote the set of all mesh points in E 2 (dots in Fig. 3(a)). Each 
(j,n) € *F is associated with a solution element, i.e. , SE (j,n). By definition, SE (j,n) is 
the interior of the space-time region bounded by a dashed curve depicted in Fig. 3(b). It 
includes a horizontal line segment, a vertical line segment, and their immediate neighbor- 
hood. 

— # 

For any ( x,t ) E SE(j>, n), u(x,<) and h(x,t), respectively, are approximated by 


u*(x,t;j,n) = u? + (ti r )?(x-Xj) + (u t )?(f-f n ) (1.7) 

and 

h*(x,t;j,n ) = (au*(x,f ;;,n), u*(x,t-,j,n)) (1.8) 

Note that (i) u”, (u x )J, and ( u t )" are constants in SE (j,n), (ii) ( Xj,t n ) are the coordinates 
of the mesh point (j, n), and (iii) Eq. (1.8) is the numerical analogue of the definition 
h = ( au , u ). 

Let u = u*(x, t ; \j , n) satisfy Eq. (1.6) within SE(j, n). Then one has (u t )” = —a (u x )J. 
As a result, Eq. (1.7) reduces to 

u*(x,t\ j,n) = u 1 - + (u x )" [(x - Xj) - a{t - <”)], (x, t) E SE(j, n) (1.9) 

i.e., u" and (u x )j are the only independent marching variables associated with (j, n). 

Let E 2 be divided into nonoverlapping rectangular regions (see Fig. 3(a)) referred to 
as conservation elements. As depicted in Figs. 3(c) and 3(d), two CEs, i.e., CE_(/,n) and 
CE+(/,n), are associated with each interior mesh point (j, n) 6 'F. These CEs will be 
referred to as basic conservation elements (BCEs). Contrarily, CE(j, n) (see Fig. 3(e)), 
which is the union of CE ~(j,n) and CE+(j,n), will be referred to as a compounded 
conservation element (CCE). 

Note that, among the line segments forming the boundary of CE _(j, n), AB and 
~AD belong to SE(j,n), while CB and CD belong to SE(j — 1/2, n — 1/2). Similarly, the 
boundary of CE-f.^’, n) belongs to either SE (j,n) or SE(j + 1/2, n — 1/2). As a result, by 
imposing two conservation conditions at each (j, n) 6 'F, i.e., 


h* - ds -- 0, 


S(C E±(j,n)) 


(j,n) € *F 


( 1 . 10 ) 
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and using Eqs. (1.8) and ( 1.9), oih' has (i) 


u " = 5 {(1 + *')V,' /j +(1 + (1 - i' 2 ) 




and, assuming 1 — i/ 2 ^ 0. (ii! 


K + )" = o 


^ - «; r , l / 2 2 - (i - - ( i + ^)( 4);; 1 1 / 2 2 


( 1 . 12 ) 


Here t/ = a&t/&x and d = f (ax/4)(ux )". The a scheme [1.5.8], the explicit nondissi- 

pative CE/SE solver for Eq. (1.6). is formed by Eqs. (1.11) and (1.12). 

According to Eq. (1.10), the total flux of h* leaving the boundary of any BCE is zero. 
Because the surface integration over any interface separating two neighboring BCEs is 
evaluated using the information from a single SE, obviously the local conservation relation 
Eq. (1.10) leads to a global flux conservation relation, i.e., the total flux of h* leaving the 
boundary of any space-time region that is the union of any combination of BCEs will also 
vanish. In particular, because CE (j,n) is the union of CE_(j,n) and CE + (j, n), 

<f h* ■ ds — 0, (j>n)eV (1.13) 

JS(CE(j,n)) 

must follow from Eq. (1.10). In fact, it can be shown that Eq. (1.13) is equivalent to 
Eq. (1.11). 

In addition to the nondissipative a scheme, there is a broad family of dissipative CE/SE 
solvers of Eq. (1.6) in which only the less stringent conservation condition Eq. (1.13) is 
assumed [2, 3, 5, 8]. Because Eq. (1.13) is equivalent to Eq. (1.11), for each of these schemes, 
u" is still evaluated using Eq. (1.11) while (u + )” is evaluated using an equation different 
from Eq. (1.12). Among these schemes is one (referred to as the a-a scheme) which is among 
the simplest and yet capable of handling solutions with discontinuities. For this scheme, 
(u+ )™ is evaluated using a finite-difference/weighted-average procedure which involves a 
parameter a (see Eqs. (2.62), (2.63) and (2.65) in [12]). The key disadvantage of the a-a 
scheme and its extensions (see below) is that, compared with the more general CE/SE 
schemes, they allow for less freedom in adjusting numerical dissipation. As explained 
in Sec. 5.5 of [9], this inflexibility may impose a constraint on the performance of these 
schemes in numerical simulations involving highly nonuniform meshes. 

The above description of the CE/SE development is based on a simple PDE. However, 
it represents the essence of the general CE/SE development which may involve a system 
of conservation laws in one, two or three spatial dimensions. In particular, note that: 

(a) The ID Euler extension of the a-a scheme, which first appears in [2], has been shown 
to be an accurate and robust shock-capturing solver [2, 3, 5. 6]. 

(b) In the original 2D extension of the CE/SE method [4,6-10], triangles are used as the 
basic building blocks of the spatial meshes. Corresponding to the three sides of a 
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triangle, three BCEs are defined for each mesh point. The union of the three BCEs at 
a rrxesh point form the CC'E at the same mesh point. Among the family of 2D CE/SE 
schemes described in [4,G 10], the 2D a scheme, which has three unknowns u, u x and 
Uy at each mesh point, are constructed by imposing three conservation conditions 
over the three BCEs at each mesh point. On the other hand, only one conservation 
condition (imposed over the CCE) per mesh point and per conservation law is used in 
the construction of the 2D Euler a-a scheme (i.e., the scheme defined by Eqs. (6.54), 
(6.107) and (6.108) in [8]). Because of its simplicity, accuracy and roubustness, all the 
numerical results presented in [4,8,9] are generated using the 2D Euler a-a scheme, 
(c) The 3D Euler a-a scheme [11] is a straightforward extension of the 2D Euler a-a 
scheme taking into account that: (i) tetrahedrons are used as the basic building 
blocks of 3D spatial meshes; and (ii) corresponding to the four sides of a tetrahedron, 
the CCE at each mesh point is the union of the four BCEs defined at the same mesh 
point. 

1.4. The Objectives and Outline of the Present Work 

In this paper, the 2D and 3D unstructured-mesh a-a Euler schemes will be constructed 
using quadrilateral and hexahedral meshes, respectively. It will be shown that the present 
schemes are also simple, robust, and accurate. The rest of the paper is organized as follows. 

The 2D and 3D solvers along with their key properties are described in Secs. 2 and 3, 
respectively. Numerical examples axe presented in Sec. 4 to demonstrate the capabilities 
of the present solvers. The concept of local and global flux conservation for the present 
2D scheme with an unstructured mesh along with a post-marching procedure for handling 
a possible “solution decoupling” problem is discussed in the Appendix. The concluding 
remarks are given in Sec. 5. 

2. The 2D Unsteady Euler Solver 

Consider the standard conservation form of the two-dimensional unsteady Euler equa- 
tions of a perfect gas [9]: 


du m dfm dfjm 

dt dx dy 


m = 1, 2, 3,4 


( 2 . 1 ) 


where f m and g m , m = 1,2, 3, 4, are explicit functions of the independent flow variables 
u m , m = 1,2, 3, 4 [9]. Let x\ = x, X 2 = y and x 3 = t be the coordinates of a three 
dimensional Euclidean space E 3 . Then, in the case that u m are smooth functions of x, y, 
z and t, Eq. (2.1) can be derived from the more fundamental conservation laws 

(p h m • ds = 0, m = 1,2, 3,4 (2.2) 

Js(V) 


where (i) S(V ) and ds were defined following Eq. (1.4); and (ii) h m (/ m ,<7 m ,u m ). Note 
that Eq. (2.2) is valid even in the presence of flow discontinuities. 
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For the; future development, let 


fm,e'= dfm/dut, g m ,f'= dg m /du e , mj. = 1,2, 3, 4 (2.3) 

2.1. Conservation Elements and Solution Elements 

Consider Fig. 4(a). Here the x-y plane is divided into nonoverlapping convex quadri- 
laterals and any two neighboring quadrilaterals share a common side. Moreover, (i) vertices 
and centroids of quadrilaterals are marked by dots and circles, respectively; (ii) Q is the 
centroid of a typical quadrilateral B l B 2 B- i B i : (iii) A x , A 2 , A3 and A 1t respectively, are the 
centroids of the four quadrilaterals neighboring to the quadrilateral B_iM_ 2 B 3 B a ; and (iv) 
Q* (marked by a cross) is the centroid of the polygon A x B_\A 2 Bj 2 A^B^ i A A B A . Hereafter, 
point Q* (which generally does not coincide with point Q ) is referred to as the solution 
point associated with the centroid Q. Note that points A*, A. 2 , A3 and A£, which are also 
marked by crosses, are the solution points associated with the centroids A 1? A 2 , A3 and 
A v respectively. 

Next consider Fig. 4(b). Here (i) t = n&.t at the nth time level (n = 0, 1/2, 1, 3/2, . . .); 
and (ii) for a given n > 0. Q. Q' , and Q" , respectively, denote the points on the nth, the 
(n — l/2)th, and the (n + l/2)th time levels with point Q (see Fig. 4(a)) being their common 
spatial projection. Other space-time mesh points, such as those depicted in Fig. 4(b), and 
also those not depicted, are defined similarly. In particular, (i) Q *, Aj, A 2 , A3 and A£, 
by definition, lie on the nth time level and, respectively, are the space-time solution mesh 
points associated with points Q , Ai, A2, A3 and A4; and (ii) Q'*, Aj*, A 2 * , A3* and A4*, by 
definition, lie on the (n — l/2)th time level and, respectively, are the space-time solution 
mesh points associated with points Q' , A\, A 2 , A3 and A4. 

With the above preliminaries, the solution element of point Q *, denoted by SE(Q*), 
is defined as the union of the five plane segments Q'Q" B"B[, Q' Q" B 2 B 2l Q' Q" B 2 B' z , 
Q' Q" B 4 B 4 , and A 1 i? ! A2 B 2 A3 By, A 4 B .\ , and their immediate neighborhoods. Moreover, 
the four basic conservation elements (BCEs) of point Q , denoted by CE^(Q), £ = 1, 2, 3, 4, 
are defined to be the space-time cylinders A\B\QB i A' l B\Q' B' A , A 2 B 2 QBiA 2 B 2 Q' B[, 
A 3 B 3 QB 2 A'^B'^Q' B 2 and A a B a QB^A\B\Q' B' 2 , respectively. In addition, the compounded 
conservation element (CCE) of point Q, denoted by CE(Q), is defined to be the space- 
time cylinder A\B\A 2 B 2 AzB^A l xBiA\B\A' 2 B 2 A' 2 B 2 A\B , i , i.e., the union of the above four 
BCEs. 

In this section, (i) the set of the space-time mesh points whose spatial projections are 
the centroids of quadrilaterals depicted in Fig. 4(a) is denoted by Q; and (ii) the set of 
the space-time mesh points whose spatial projections are the solution points depicted in 
Fig. 4(a) is denoted by fi*. Note that the BCEs and the CCE of any mesh point € and 
the SE of any mesh point 6 Q* are defined in a manner identical to that described earlier 
for point Q and Q*. 

2.2. Approximations Within a Solution Element 

For any Q* € ft* and any (x,yj.) G SE(Q*), u m (x,y,t), ffm(x,y,t) and 

h m (x, y, t), respectively, are approximated by u* m {x,yA\ Q *), fm( x , ?/A; Q *), £^(x, y,t; Q*) 


and h^(x, y, t\ Q*) (see below). For any m ~ 1,2, 3, 4, let 

u* m (x.yJ.;Q*) = f (u m ) Q . + (u mx ) Q ‘{x - x Q . ) + (u my ) Q . (y - ijq. ) + (u ml ) Q .(t - t n ) (2.4) 

where (i) (xQ«,yg-,f n ) are the coordinates of the space-time solution mesh point Q *; 
and (ii) (u m )g-, ( it m :t)g ., ( u mj ,)g * and {u ml )Q-, which are constants in SE(Q*), are the 
numerical analogues of the values of u m , du m /dx. du m /dy and du m /di at point Q * , 
respectively. 

Let (fm)Q-, (gm)Q-, ( fm,e)Q * and ( g m ,e)Q • denote the values of the functions / m , g m , 
f m ,e and g m j , respectively, when u m , m = 1,2, 3, 4, respectively, assumes the values of 
(a m )Q*, m = 1,2, 3, 4. Then, for any m, we define 


(Jmi)Q' = (jmi)g- = ^(gm,* )q« ( Ujx )q* 

e=i 


e=i 

4 


(fmy )q* ~ ( u ly)Q * i ( 9my)Q m ~ * ( u (y)Q m 

e=i e=i 

4 4 

(fmt)Q* d = (9mt)Q* *= {u£t)Q- 


e=i 


£= 1 


Because (i) 


3/n 


= E f*» 


dut 


dx f - dx 

e=i 


(2.5a) 

(2.56) 

(2.5c) 

(2.6) 


and (ii) the expression on the right side of the first equation in Eq. (2.5a) is the numerical 
analogue of that on the right side of Eq. (2.6) at point Q*, (/ mx )Q* can be considered as 
the numerical analogue of the value of df m /dx at point Q* . Similarly, (g mx )q* , (fmy )Q - , 
( 9my)Q * , (/mt)Q* and (g m t)Q- can be considered as the numerical analogues of the values 
of dg m /dx, df m /dy , dg m /dy , df m /dt and dg m /dt at point Q*, respectively. As a result, 
for any m = 1, 2, 3, 4, we define 


def 


/*(x,y,t;Q*) = U m ) Q . + {f mx ) Q .( x -x Q .) + {f my ) Q .(y-y Q .) + (f mt ) Q .{t-t n ) (2.7) 


and 

def 

g* m (x,y,t\Q*) = { 9 m)Q- +(gmx)Q-(x -XQ-) + (g my ) Q *(y -y Q -) + (g mt ) Q -(t-t n ) (2.8) 
Also, as an analogue to h m = (/ m , gm,u m )< for any m = 1,2, 3, 4. we define 


h* m (x,y,t\Q*) d = (fn(x,y,t\Q*),g* m (x,y,t;Q*),u* m (x,y,t\Q*)) 


(2.9) 
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Note that, by their definitions: (i) (/m)g*, (<Jm )g- , (f m> i)Q m , and (y m ,e)Q- are functions of 
(a m )g-,m = 1 , 2 , 3. 4; (ii)(/ m x)g- and (y mi )g. are functions of (u m )g- and(u mx )g., m = 
1,2,3. 4; (iii) (f my )Q- and (g my )Q' are functions of (u m )Q. and (u my )g., m = 1,2, 3, 4; 
and (iv) (/ m «)g- and (y m( )g- are functions of (u m )g. and (tt m< )g., m = 1,2, 3,4. 

To proceed, we also assume that, for any (x,y,t) £ SE(Q*), and any m — 1,2, 3,4, 


du* m (x,y,t\Q*) + df^(x,y,t;Q m ) + dg^j^yj; Q*) _ Q 


<9* 


<9x 


dy 


( 2 . 10 ) 


Note that Eq. (2.10) is the numerical analogue of Eq. ( 2 . 1 ). With the aid of Eqs. (2.4), 
(2.7), (2.8), (2.5a). and (2.5b). Eq. (2.10) implies that, for any m = 1,2, 3, 4, 


(a m e)g* — (/milo" ( 9my)Q’ — ^ ^ [(/m,l)g* ( u (x)q * 4" (dm, l )Q' { u fy )g* ] (2-11) 

e=i 

Thus (u m ()g- is a function of (it m )g», (u mi )g-, and (u mj ,)g-, m = 1,2, 3, 4. From this 
result and the facts stated following Eq. (2.9), one concludes that the only independent 
discrete solution variables associated with the space-time solution point Q* are (u m )g*, 
(«mi)Q*, and (limy )g* , m = 1 , 2 , 3, 4. 

2.3. Evaluation of (u m )g* 

Based on Figs. 4(a) and 4(b), we introduce the following preliminaries: 

(a) The boundary of CE(<Q) belongs to the union of SE(Q*), and SE(T^*), £ = 1 , 2 , 3, 4. 
Specifically, (i) the octagon A 1 B 1 A 2 B 2 A 3 B 3 A 4 B 4 belongs to SE($*); (ii) the quadri- 
laterals A\B[Q' B' 4 , A\B[B 4 Ai, and A\B[B\A\ belong to SE(A , 1 *); (iii) the quadri- 
laterals A 2 B 2 Q' B[, A' 2 B[B x A 2 , and A 2 B 2 B 2 A 2 belong to SE^^*); (iv) the quadrilat- 
erals A'^B-iQ' B' 2 , A 3 B 2 B 2 A 3 , and A 3 B 2 B 3 A 3 belong to SE(A 3 *); and (v) the quadri- 


(b) 


laterals A' 4 B' 4 Q' B' 3 , A 4 B 3 B 3 A 4 , 


and A 4 B' 4 B 4 A 4 belong to SE(^. 4 *). Note that, by 


definition, (i) the quadrilaterals A\B\QB 4 , A 2 B 2 QB 1 , A 3 B 3 QB 2 , and A 4 B 4 QB 3 , 
which form the octagon A\B\ A 2 B 2 A 3 B 3 A 4 B 4 (the top face of CE(Q)), also be- 
long to SE(T*), SE(A 2 ), SE(A 3 ), and SE(T 4 ), respectively; and (ii) the octagon 
A\ B[A' 2 B' 2 A 3 B' z A\ B ' 4 (the bottom face of the CE(Q)) also belongs to SE(Q'*). How- 
ever, in the evaluation of Eq. (2.13) (see below), by assumption, the top face of CE(Q) 
is considered to be a subset of SE(Q*) while the bottom face is considered to be the 
union of subsets of SE(A' e *), £ = 1 , 2 , 3, 4. 

Let r be a space-time plane segment lying within SE(Q*). Let (i) A be the area 
of T; (ii) ( x c ,y c ,t c ) be the coordinates of the centroid of T; and (iii) n be a unit 
vector normal to T. Then, because u* m (x,y,t-,Q*), f^(x,y,t;Q*) and y^(x, y, t; Q*) 
are linear in x, y and t, Eq. (2.9) implies that 


l 


^m( ^ a Vci^C) Q ) ' An 


( 2 . 12 ) 


where ds = do n with do being the area of a surface element on T . 
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(c) Let, S denote the area of the top face Ai D\ A 2 B >A 3 B:\A.\B.\ of CE (Q). Because 
the unit outward normal vector (outward from the interior of C'E( Q)) of this face is 
(0, 0, 1), its surface vector (i.e., the unit outward normal vector multiplied by the area) 
is (0,0,5). 

(d) Let (i) (x f ,y e ), £ = 1,2, 3, 4, denote the spatial coordinates of the centroids of the 
quadrilaterals A\B\Q' B\, A' 2 B 2 Q' B [ , A 3 B' 3 Q'B 2 and A' 4 B' 4 Q' B' 3 , respectively; and 
(ii) S e , £ = 1,2, 3. 4, denote the areas of the above four quadrilaterals, respectively. 
Then (i) (x e , y e , t n ~ 1 / 2 ), £=1,2 ,3,4, are the coordinates of the above four centroids, 
respectively; and (ii) (0, 0, —S ( ), £=1,2, 3, 4, are the surface vectors of the above four 
quadrilaterals, respectively. Furthermore, because (i) the above four quadrilaterals 
form the bottom face of CE(Q): and (ii) the area of the top face of CE(Q) is identical 
to that of the bottom face, one concludes that 5 = J2e=i 

(e) Let the eight side faces A[B' 4 B 4 Ai, A\B[B X A X , A' 2 B[BiA 2 , A! 2 B 2 B 2 A 2 , A' 3 B' 2 B 2 A 3 , 
A' 3 B' 3 B^Az, A 4 B' 3 B 3 A 4 and A' 4 B' 4 B 4 A 4 of CE(Q) be assigned the indices (1, 1), (2, 1), 
(1,2), (2,2), (1,3), (2,3), (1,4) and (2,4), respectively. Hereafter each side face with 
the indices (&,£) is referred to as the (fc,£) side face. For each £, by definition, the 
(1,£) and (2,£) side faces belong to SE(i4^*). Because (i) the spatial projection of 
each side face is a line segment on the x-y plane; and (ii) each side face is sandwiched 
between the (n — l/2)th and the nth time levels, one concludes that, for the (fc,£) 
side face, its surface vector and the coordinates of its centroid, respectively, are given 
by (At/2)\ e k (n e kx ,n e ky ,0) and (x[,y e k ,t n - a£/4). Here A£, ( n e kx ,n e ky ), and ( x e k ,y e k ), 
respectively, denote the length, the unit outward normal (on the x-y plane), and the 
coordinates of the midpoint of the spatial projection of the (k, £) side face. 

(f) Note that: (i) (xQ.,yQ-,i n ) are the coordinates of the centroid Q* of the top face 
A x BiA 2 B 2 AzBzA 4 B a of CE(Q); (ii) u* m {xQ . , y Q - , t n \ Q*) = (u m )g- (see Eq. (2.4)); 
and (iii) the surface vector of the top face is (0, 0, 5). As a result, Eq. (2.9) and (2.12) 
imply that the flux of h * m leaving CE (Q) through its top face is (u m )g* 5. Similarly, by 
using the information presented in items (a), (b), (d) and (e), the flux of h ^ leaving 
the other faces of CE(Q) can be evaluated in terms of the independent marching 
variables at points A' ( *, £ = 1,2, 3, 4. 

Let 

( l h* m -ds= 0, m = 1,2, 3,4 (2.13) 

JS{CE(Q)) 

i.e., the total flux of h* m leaving CE(Q) through its boundary vanishes. Then, with the aid 
of the above preliminaries, it can be shown that 


(u m )«■ = (Y, R l )/S. m = 1,2. 3,4 (2.14) 

e=i 


10 


where, for any ra, Z = 1,2, 3, 4, 


R f m =S , u* m (x t ,y',t n - l ' 2 ;A'n 

~ ^ IT^k i n kx fm&k'Vk't 71 ~ A */4; A ' t * ) + g* m {x k ,y e k ,t n — Af/4; A'/)] 


*=1 


Because, by definition, t = t n ~ 1 / 2 for any point A'f, here the functions u* m (x, y,t; A' e *), 
/^,(x, y, t; .4' £ *), and g* m (x, y,t\ A'f) are defined using Eqs. (2.4), (2.7), and (2.8), respec- 
tively, with the symbols Q* and t n in these equations being replaced by A and < n-1 / 2 , 
respectively. As a result, each R e m and therefore each (u rn )Q-, an independent marching 
variable at the nth time level, is a function of several independent marching variables at 
the (n — l/ 2 )th time level, i.e., (u m )/i'* , ( Umx)A '* and i u my)/ t'* , rn,£ = 1 , 2 , 3, 4. 

2.4. Evaluation of (u m x)<3’ and {u my )Q- 

A finite-difference approach similar to that given in [10] is employed here to evalu- 
ate (u mI )Q- and ( u my )Q *. First, we perform a spatial translation of the quadrilateral 
A* A% A 3 A^ so that the centroid of the resulting new quadrilateral A^A^A^A^ coincides 
with Q* (see Fig. 4(c)). Let the centroid of the quadrilateral A*A^A^A^ and its spatial 
coordinates be denoted by A * and (x a- , y a*), respectively. Then (xa the spatial 

coordinates of A\, are 


xa° = x^; + xq* - xa-, and y A ° = Va, ; + VQ- ~ JM* , £=1,2, 3,4 

To proceed, let 




m,l = 1 , 2 , 3, 4 


(2.16) 


(2.17) 


Next, for any m = 1,2, 3, 4, consider the three points in the x-y-u space with the coor- 
dinates (xq* , j/q* , (u m )Q- ), (x A o,y j 4 »,(u m ) /1 ») and (x A ° 2 , yA° 2 , («m )/!•), respectively. The 
values of du/dx and du/dy on the plane that intercepts the three points are given by 


where 


and 


(«21)«-“a,/A and (u^)q- d = A y /A (A + 0) (2.18) 


A x 


A = f 


*A° 


def 


( )/t° 
{ u m) A° 


- xq- Va\ 

- XQ- Va\ 

( u m)Q * 
( u m )Q m 


~ VQ * 

~ VQ * 

yA\ - yQ- 
yA° 2 - yQ- 


(2.19 a) 
(2.19 b) 


A 


y 


def 


( u m) A° (Wjji )q* 

(a m) A 2 ( u m )q* 


xq* — X^o 

XQ- — X.4» 


(2.19c) 
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Note that: (i) A = 0 if and only if the spatial projections of .4", A% and Q* are collinear; 

and (ii) similarly, (uml)g* and (u^)q. , k = 2,3,4, are defined, respectively, by replacing 
the points A" and A£ in the above operations with (i) A£ and A£; (ii) and A 4 ; and (iii) 
A^ and A°, respectively. 

With the above preliminaries, for each m = 1,2, 3, 4, (u mx )g* and (u mj ,)g. may be 
evaluated by 


1 4 1 4 

( u mx)g* = “ )<?* » ( u mj/)g* = ~ (2.20) 

k— 1 A:=l 

Alternatively, for a flow with steep gradients or discontinuities, the simple averages in 
Eq. (2.20) may be replaced by weighted averages, i.e., 


(^mi)g* — 


1°’ 


1 EL, 

(Wi k) ) a (u ( ^l) Q .' 


if Qmk = 0,k = 1,2, 3, 4 


and 


( U my)g* — 


f°, 


i EL , 



if &mk = 0, k = 1,2, 3, 4 


(2.21a) 


(2.216) 


Here (i) a > 0 is an adjustable constant (usually a = 1 or a = 2); (ii) 


a def 
v mk — 


(Umi)g. + (timy)g* , m,k = 1,2, 3, 4 


( 2 . 22 ) 


and (iii) 

= f o m2 6 m3 e m4 , e m3 e m4 e mU w. i 3) d = e m4 9 ml 9 m2 , w. i 4) d = e ml $ m2 0 m 3 

(2.23) 

Note that: (i) to avoid dividing by zero, in practice a small positive number such as 10 -60 
is added to the denominators that appear in Eqs. (2.21a) and (2.21b); and (ii) Eqs. (2.21a) 
and (2.21b) reduce to Eq. (2.20) if a = 0. 

2.5. Remarks and Discussions 

The present 2D Euler solver is formed using Eqs. (2.14), (2.21a) and (2.21b). Stability 
of the solver generally requires that (i) a > 0, and (ii) the maximal CFL number < 1. 
Also, (i) with a > 1, the solver is capable of suppressing numerical oscillations near a 
discontinuity; and (ii) solutions generated by the solver tend to become more smeared as 
the CFL number decreases or the value of a increases. Other key properties of this solver 
are given in the following remarks: 
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(a) The stencil of the present explicit, solver is formed by one point at the upper time level 
and four points at the lower time level. Because the spatial projections of the four 
points at the lower time level are the immediate neighbors of that of the point at the 
upper time level, the stencil is staggered in space-time, and it is most compact among 
the schemes using quadrilateral meshes. As a result, the solver is ideal for parallel 
computations. 

(b) For a uniform mesh, points such as Q. Q* and .4* referred to earlier coincide with one 
another. In this case, the present solver can be greatly simplified. Also, by using the 
arguments presented in [8,9] and also by numerical experiments, it can be shown that 
the simplified scheme is second order in accuracy. 

(c) The present scheme is applicable to both structured and unstructured meshes. For 
a structured mesh, the set Q* may be divided into two disjoint subsets and 
with the following property: If any point, say point Q* , belongs to fi!j_ (fil), then the 
six space-time solution mesh points immediately neighboring to point Q*, i.e., points 
Q '* , Q"*, and A* e , £ = 1,2, 3,4, belong to Cl*_ (Q+). Because, for each £ = 1,2, 3,4, 
points A' ( * and A* t are immediate neighbors of each other and thus they must belong to 
different subsets, one concludes that points Q*, and A' e *, L = 1,2, 3, 4, which form the 
stencil of the present marching scheme, belong to the same subset. From the above 
observations, it is seen that each of and Q.*_ represents a staggered space-time 
mesh. As such, the entire space-time mesh is a dual space-time mesh [9], i.e., the 
union of two disjoint staggered space-time meshes. Furthermore, it is also obvious 
that the marching over is completely decoupled from that over Q*_, i.e., marching 
needs to be carried out only over one of these two staggered space-time meshes, unless 
the decoupling is prevented by other factors such as the boundary conditions imposed. 
Note that boundary values generally are not updated using the main marching scheme. 
As a result, solution values of and Q*_ may become coupled near a boundary (see 
Sec. 4). 

(d) Consider the decoupling case referred to in item (c). Let a space-time mesh point 

belong to Q+ (fi_) if and only if its associated space-time solution mesh point belongs 
to (QI). Then it is obvious that the set Q is formed by the two disjoint sets 

and Q_. Moreover, the CCEs of the mesh points in (f)_) do not overlap 

among themselves and they can fill any domain in E 3 . Furthermore, because the 
surface integration over any interface separating two neighboring and nonoverlapping 
CCEs is evaluated using the information from the same SE (i.e., the flux leaving a 
CCE through its interface with a neighboring CCE is the negative of the flux leaving 
the neighboring CCE through this interface), a summation of the local conservation 
conditions Eq. (2.13) over the mesh points Q € Q+ (£2_) leads to a global conservation 
condition, i.e., for each m = 1,2, 3,4, the total flux of h* m leaving the boundary of 

any space-time region that is the union of any combination of the CCEs associated 

with fl + (fi_) vanishes. Note that a similar discussion for the general case in which 
decoupling may not occur will be given in the Appendix. 

(e) The present solver and the triangular-mesh-based solver described in [10] are con- 
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structed using similar techniques. Using these techniques and their trivial extensions, 
one can easily develop a 2 D CE/SE solver for spatial meshes formed by polygons of 
different shapes. An advantage of using such a mixed mesh is that a geometrically 
complex spatial subdomain can be filled easily using triangles while a less complex sub- 
domain, such as a near-wall region, can be filled using more regular shaped polygons 
such as quadrilaterals. 

(f) Because of the space-time staggering nature of the stencil of the present scheme, a 
solution of the present scheme may appear as the overlapping of two distinctively dif- 
ferent solutions (especially in a high-gradient region) after many marching steps. The 
significance of this Solution decoupling'' problem and how to handle it are discussed 
in the Appendix. Note that this problem could occur even in the absence of a complete 
decoupling referred to earlier Also because the solution decoupling problem 
is not significant for the test problems discussed in Sec. 4, the numerical results pre- 
sented there are generated without using the post-marching procedure described in 
the Appendix. 

3. The 3D Unsteady Euler Solver 

For the current 3D case, Eqs. (2.1 )— (2.3) are replaced by 


and 


^Urn dfm 9g m 

dt dx dy dz 


m = 1 , 2 , 3 , 4, 5 


(p h m • ds = 0 , m = 1 , 2 , 3, 4, 5 

Js{V) 


(3.1) 

(3.2) 


fm,t = dfm/dut, g m , t = dg m /du t , qm,i = dq m /du e m,£ = 1,2, 3, 4, 5 (3.3) 

-* ^ef 

respectively. Here (i) h m — i and (ii) the three dimensional Euclidean 

space E 3 referred to in Sec. 2 is replaced in the current case by the four dimensional 
Euclidean space E 4 with x\ = x, X 2 = y , £3 = z , and £4 = t. 

3.1. Conservation Elements and Solution Elements 

The spatial computational domain is divided into nonoverlapping convex hexahedrons 
of arbitrary shape with the understanding that any two neighboring hexahedrons share a 
common face. In Fig. 5, Q (marked by a circle) is the centroid of a typical hexahedron 

(hereafter referred to as the central hexahedron). Each of the 
central hexahedron's six neighboring hexahedrons is arbitrarily assigned an identification 
index £ = 1 , 2 ,..., 6, i.e., the neighboring hexahedron with the index £ is referred to as 
the ^th neighbor of the central hexahedron. Also the centroid of the ^th neighbor will be 
denoted by A As an example, the central hexahedron and its first neighbor is separated 
by the quadrilateral in Fig. 5. 

With the above preliminaries, we proceed with the following definitions: 
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(a) Point Q and the two end points (say points B_\ and Z? 2 ) of any of the twelve edges 
of the central hexahedron form a triangle. Each of the twelve triangles so formed is 
arbitrarily assigned an index j = 1, 2, 3, . . . , 12 and denoted by A (j). 

(b) Given any Z = 1,2, . . . , 6, a triangle is formed by the point and the two end points 
of any of the four edges of the interface (a quadrilateral) that separates the central 
hexahedron and its Zth. neighbor. Each of the four triangles so formed with the same Z 
is arbitrarily assigned an index k = 1.2, 3, 4. and denoted by A(k,Z). As an example, 

Ad^,^, A and A A l B_ 5 B_ 1 depicted in Fig. 5 have the same 
Z = 1. Therefore they may be denoted by A(l,l), A(2, 1), A(3, 1), and A(4, 1), 
respectively. 

(c) The centroid of the 24-faced polyhedron B_[B_ 2 ^ 3 MaS.^B_ 6 B_ 7 B s AiA 2 A 3 A 4 A 5 A 6 is 
referred to as the solution point associated with point Q. Note that (i) the above 
24- faced polyhedron hereafter is denoted by F(24); and (ii) the centroid of V/24) is 
denoted by Q* and marked by a cross in Fig. 5. 

(d) Given any Z = 1,2,..., 6, points Q , A^ and the four vertices of the quadrilateral 
interface that separates the central hexahedron and its Zt h neighbor are the vertices 
of a octahedron. This octahedron hereafter is denoted by V_(8;Z). 


In the space-time computational domain, again we assume that t = n&t at the nth 
time level (n = 0, 1/2, 1,3/2, . . .). Also, for a given n > 0, let Q , Q ' , and Q" (not shown), 
respectively, be the points on the nth, (n — l/2)th, and (n + l/2)th time levels with point 
Q being their common spatial projection. Other space-time mesh points such as (i) Q* and 
Q'*; (ii) B k , B' k , and B%, k = 1, 2, 3, 4, 5, 6, 7, 8; and (iii) A t , A* ( , A' t and A' e *,Z= 1,2,..., 6, 
are defined similarly. Because geometric objects in E 4 generally are difficult to visualize, 
they will be described analytically in the following discussions. 

To proceed, note that a “plane” (termed a hyperplane) in E 4 , by definition, is a 
subspace of E 4 defined by a linear equation, i.e., 


cqx + a 2 y + a z z + a 4 t = a 0 ((a!) 2 + (a 2 ) 2 + (a 3 ) 2 + (a 4 ) 2 ^ 0) (3.4) 


where a*, k = 0,1,2, 3,4, are constants. As a result, a hyperplane in E 4 is a three 
dimensional subspace. The unit normal to the hyperplane is 


- _ _|_ (eg, 02 , 03 , 04 ) / 3 ^ 

V(ai ) 2 + (a 2 ) 2 + (a 3 ) 2 + (a 4 ) 2 

Note that a hyperplane segment, by definition, is a bounded region of a hyperplane. 

Two types of hyperplane segments in E 4 are involved in the definition of SEs to be 
given shortly. A hyperplane segement of type I, denoted by r(F;f c ), is formed by all the 
points (x,j/,z,t) that satisfy the conditions (i) t = t c ; and (ii) (x,j/,z) 6 V, where t c is 
a constant and V_ denotes a 3D spatial region. Obviously the equation t = t c is a special 
form of Eq. (3.4). Also it can be shown that: 

(a) The unit normal to T(F;t c ) is (0, 0, 0,±1). 

(b) The “area” of T(F;t c ) is the volume of V. 
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(<■-) The coordinates of the centroid of T( VjJ r ) are (x c , y„, z c , t r ) where (x r ,y r ,z r ) are the 
coordinates of the centroid of V. 

On the other hand, a hyperplane segment of type II, denoted by T(5; <_,<+), is formed 
by all the points (x, y , 2 , t) that satisfy the conditions: (i) (x, y, z) G 5; and (ii) f_ < t < < + , 
where 5 denotes a spatial plane segment, and <_ and t + (f_ < t + ) are constants. Note 
that every point (x,y,z-) on the spatial plane segment 5 satisfies a linear equation of the 
form 

c\x + c 2 y + c 3 z = c 0 ((c x ) 2 + (c 2 ) 2 + (c 3 ) 2 ^ 0) (3.6) 

where c*, k = 0,1, 2, 3, are constants. Thus every point ( x,y,z,t ) on r(5;t_,t + ) also 
satisfies a special form of Eq. (3.4), i.e.. Eq. (3.6). Moreover, it can be shown that: 

(a) The unit normal to r(S;f_,f-t-) is (n,0) where n is the unit normal to the spatial 
plane segment 5, i.e., 

- _ _j_ (ci,c 2 ,c 3 ) 

a/( c i ) 2 + (c 2 ) 2 + (c 3 ) 2 

(b) The “area” of r(5;f_,f+) is the area of 5 multiplied by ( t+ — f_). 

(c) The coordinates of the centroid of r(5;£_,t+) are (x c , y c , z c , (t_ + t + )/2) where 
(x c , y c , z c ) are the coordinates of the centroid of 5. 

In addition to the above two types of hyperplanes, we shall also consider “hyper- 
cylinders” in E\. A hypercylinder, denoted by A (F;t_,t + ), is formed by all the points 
( x,y,z,t ) that satisfy the conditions: (i) (x, y, z) 6 V; and (ii) <_ < t < t + , where F is a 
3D spatial region, and t_ and (f_ < t+) are constants. 

With the above preliminaries, SE(<5*), the solution element of point Q * — the point 
that lies on the nth time level and has Q* as its spatial projection, is defined to be the 
union of T(F(24); f n ) and T(A(y); f n-1 / 2 , < n+1 / 2 ), j = 1, 2, 3, . . . , 12, and their immediate 
neighborhoods. Moreover, the six basic conservation elements (BCEs) of point Q , denoted 
by CE*(Q), £ = 1,2, ..., 6 , are defined to be the hypercylinders A(F(8;£);f n-1 / 2 ,t n ), 
i = 1,2, ...,6, respectively. In addition, the compounded conservation element (CCE) of 
point Q, denoted by CE(<5), is defined to be A(F(24); f n-1 / 2 , t n ), i.e., the union of the 
above six BCEs. 

In this section, (i) the set of the space-time mesh points whose spatial projections are 
the centroids of the hexahedrons that fill the 3D spatial computational domain is denoted 
by fi; and (ii) the set of the space-time mesh points whose spatial projections are the 
solution points of the centroids referred to in item (i) is denoted by Q,*. Note that the 
BCEs and the CCE of any mesh point G and the SE of any mesh point € are defined 
in a manner identical to that described earlier for point Q and Q*. 

3.2. Approximations Within a Solution Element 

For any Q* € ft* and (x,y,z,t) G SE(Q*), u m (x,y,z,t), f m (x,y,z,t), y m (x,y, 2 ,t), 
q m (x,y,z,t), and h m (x,y,zj) are approximated by u* m (x,y,z,t-,Q*), f^(x,y,z,t]Q*), 
g* m (x , y, z,t; Q*), q„(x, y,z,t: Q*), and h^(x,y, z,t;Q*), respectively (see below). For any 


(3.7) 
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m = 1, 2, 3, 4, 5, let 


def 


i(x, y. z, t\ Q*) = (u m ) Q . 4- («mi)Q'(l — XQ- ) 4- (Umy)Q'(y ~ VQ' 
4- (Umz)Q*(z ~ ZQ- ) + ( U ml )q-(t — <") 


(3-8) 


f? n (x, y,Z,t;Q*) = (/m)Q* + {fmz)Q‘(x — Xq* ) + (/mj/)g*(y — UQ') 
+ {fmz )q‘ ( Z — Zq- ) 4 - (/mt)< 3 *(^ — t ) 

def 

Q ) = (,9m)<3* “1“ (Smx )q* ^Q* ) “b (<7m y )q* ( J/ J/Q* ) 

“1“ (Smz)Q*(* )"b($ f m*)g*(^ t ) 

def 

) = (?m)g* +(9mx)Q*(^ ~ ^Q* ) + (<?my )q* (j/ “ J/Q* ) 

+ — Zq * ) + (<Zm*)<2*(* — ^ ) 

and 


(3.9) 

(3.10) 

(3.11) 


h' m {x, y, 2, f; Q') d = J, 2, t; <?'), 9^(1, !/, 2, t; O’), jr, 2, f; < 3 *), <,(2, t/, 2, f; Q")) 

( 3 - 12 ) 

be the 3D extension of Eqs. (2.4) and (2.7)-(2.9). Note that, in this section it is implicitly 
assumed that any notation that has a similar 2D version is defined similarly. The definition 
of such a notation will not be given explicitly here unless confusion could occur. 

Moreover, we assume that, for any ( x,y,z,t ) 6 SE(Q*), and any m = 1,2, 3, 4, 5, 


du* m (x, y, t; Q^)_ + df^(x,y,t\Q*) + dg^(x,y,t;Q”) + dq„(x, y, t; Q m ) 


dt 


dx 


dy 


dz 


= 0 (3.13) 


Thus, for any m = 1, 2, 3, 4, 5, 


( u mt)Q * — {fmx )q* ( Qmz)Q m 

5 

= - [(/m,/)g* (utx)Q- + ( 9m,e)Q * («<») Q* + (^m.^Q* ( u Iz)q' ] 

t- 1 


(3.14) 


Using the equations given above, it can be shown that, for the current 3D case, the only 
independent discrete variables associated with the space-time solution point Q* are (u m )g* , 
(^mi)(J* ) {^my)Q m i &nd (v> mz )Q* , 171 = 1 , 2 , 3 , 4 , 5 . 

3.3. Evaluation of (u m )Q* 

We begin with the following preliminaries: 

(a) The boundary of CE(Q) is formed by the “top face” r(F(24);£ n ), the “bottom face” 
r(V(24); t n-1/2 ), and the 24 “side faces” T( A(Jfc, £)\ T" 1 / 2 , t n ), k = 1 , 2 , 3, 4 and £ = 
1,2,..., 6. Because V(24) is the union of V_(8;£), £ = 1,2, ...,6, the top (bottom) 
face is the union of r(F(8; £)\ t n ) (r(F(8; £): t n ~ 1 / 2 )), £ = 1, 2, . . . , 6. From the above 
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observations, one concludes that the boundary of C'E( Q ) belongs to the union of 
SE(Q*), and SE(A'*), 1 = 1,2,..., 6. Specifically, (i) T(V(24):t n ) belongs to SE (£•); 
and (ii) for each £ = 1,2,. .., 6 , T( F(8; £); f" -1 / 2 ) and T(A (k,£);t. n ~ l / 2 ,t n ), k = 
1,2, 3. 4, belong to SE(A'/). Note that: (i) T(V/8; £)\ t n ), £ = 1,2,..., 6, the union 
of which is T(J7(24); t n ), also belong to SE(A^), £ = 1 , 2, .... 6, respectively; and (ii) 
r(H(24); / n - |/2 ), which is the union of T(F(8; £); t n ~ ’/ 2 ), £ = 1, 2, . . . , 6, also belongs 
to SE(Q'*). However, in the evaluation of Eq. (3.16) (see below), by assumption, 
r(V( 24 ); t n ) is considered to be a subset of SE( Q * ) while T(F(24); t n ~ 1 ^ 2 ) is considered 
to be the union of subsets of SE(A' ( *), £ = 1, 2, . . . , 6. 

(b) Let T be a hyperplane segment lying within SE(Q*). Let (i) A be the area of T; 
(ii) (x c ,y C i z c ,t c ) be the coordinates of the centroid of T: and (iii) n be a unit vector 
normal to F . Then it can be shown that 

J h* m • ds = h* m (x c ,y c ,z c ,t c ;Q*) • An (3.15) 

where ds = da n with da being the area of a surface element on T. 

(c) Let V denote the volume of V_(2A), i.e., the area of the top face T(F(24); t n ) of CE(Q). 
(see comments (a)-(c) given following Eq. (3.5)). Because the unit outward normal 
vector (outward from the interior of CE((J)) of this face is (0, 0, 0, 1), its surface vector 
(i.e., the unit outward normal vector multiplied by the area) is (0,0,0, V). 

(d) Let V e and (x t ,y t ,z l ), repectively, denote the volume and the spatial coordinates 
of the centroid of any 1L(8;£). Then the surface vector, and the coordinates of the 
centroid of T(y(8; £); t n ~ 1 ^ 2 ), respectively, are (0,0,0,— V e ) and (x e ,y e , z e 

(e) Let S e k , (n[xi n iyi n kz)i an( ^ ( x ii Vk' z k)' respectively, denote the area, the spatial unit 
outward normal, and the coordinates of the centroid of any A (k,£). Then the sur- 
face vector, and the coordinates of the centroid of the side face T(A(k, 
respectively, are (At/2)S ( k (n kx , n e ky , n kz , 0) and (x k ,y k ,z k ,t n — At/4) (see comments 
(a)-(c) given following Eq. (3.6)). 

(f) Note that: (i) (xq- ,yq- , zq» ,t n ) are the coordinates of the centroid Q* of the top 

face r(V(24),t") of CE(Q); (ii) u* m (x Q - ,y Q - , zq. ,t n \Q*) = (u m )Q* (see Eq. (3.8)); 

and (iii) the surface vector of the top face is (0,0,0, V). As a result, Eq. (3.12) and 

(3.15) imply that the flux of h ^ leaving CE(Q) through its top face is (u m )Q-V r . 

Similarly, by using the information presented in items (a), (b), (d) and (e), the flux 
— # 

of h * m leaving the other faces of CE(Q) can be evaluated in terms of the independent 
marching variables at points A'f , £ = 1,2, 3, 4, 5, 6. 

Let 

<p h* m • ds = 0, m = 1, 2, 3, 4, 5 (3.16) 

JS(CE(Q)) 

i.e., the total flux of h ^ leaving CE(Q) through its boundary vanishes. Then, with the aid 

of the above preliminaries, it can be shown that 

6 

(«m)Q. =(J2 R m)/ V ’ rn = 1,2, 3,4,5 (3.17) 

e=\ 
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where, for any m = 1, 2, 3, 4, 5 and any £ = 1, 2, . . . , 6, 




pf \/t * / J X £ in— 1 / 2 . a i* \ _ V~^ q£ „£ r* / i n \ + /A- A h 

/t m v i y i z it i ^ c\ ^k ^kx /m' 'X' k i !Jk ’ ^ k * ^ A^/4, 


Jt=l 


(3.18) 


+ n i» 9m( x ki Vk - — a V4: -4^*) + yi- ~i'^ n ~ 4; -4^*) 


Here u„(x,y, z,t; A' ( *), f^(x,y, z,t; A' ( *), g^(x,y, z,t; A' e *), and q^{x,y, z,t; A' e *) are de- 
fined using Eqs. (3.8)— (3.11), respectively, with the understanding that the symbols Q* 
and t n in these equations be replaced by A£* and t n-1 / 2 , respectively. As a result, each 
R e m and therefore each (u m )q. , an independent marching variable at the nth time level, 
is a function of several independent marehing variables at the (n — l/2)th time level, i.e., 
( Um)A \ *, (u my )A' t m , and (u mz ) A >., m = 1,2, 3, 4, 5 and t = 1,2, ..., 6 . 

3,4. Evaluation of ( u mx )Q - •> («mj,)Q* and ( u mz )Q - 

First, we perform a spatial translation of the polyhedron .A* Ag Ag so that 

the centroid of the resulting new polyhedron A° A^AgA^AgAj? coincides with Q* . Let (i) 
the centroid of the polyhedron A’AjAgA^AgAg and its spatial coordinates be denoted 
by A* and (x A * , y A * , z A * ), respectively; and (ii) Sx = Xq> — x A -, Sy = yQ - — y A », and 
6z = zq • — z A - . Then (x A °,y A °,z A o), the spatial coordinates of A° t , t = 1,2, ...,6 are 
given by 


x A o - x A - ( + Sx, y A o = y A - + Sy and z A o = z A * + 6z (3.19) 

As a preliminary for the following discussions, for m = 1,2, 3, 4, 5 and i = 1, 2, . . . , 6, let 


= u* m {x A o,y A o,z A o,t n \ A' e *) 


(3.20) 


Su 1 ^ ( Um ) A o ~{u m ) Q . (3.21) 

and 

6 xt = x A o -xq., Sy e = y A o -y Q *, Sz t = z A o - z Q * (3.22) 

Next consider the vertex B_ x depicted in Fig. 5. This vertex is the common vertex of 
the central hexahedron and three of its neighbors. As an example, let the identification 
indices £ of these three neighbors be 1, 2, and 3. Then, for any m = 1,2, 3, 4, 5, consider 
the four points in the x-y-z-u space with the coordinates ( xq * , yQ* , zq* , (um)Q m ), and 
(x A °,y A °,z A o,(u m ) A o), i = 1,2,3. It can be shown that the values of du/dx, du/dy , and 
du/dz on the hyperplane that intercepts the above four points are given by 


(u 


( 1 ) 

TUX 


\ def a 

)q. = A* 


/A, 


where 



A d = 


Ay/A, 

(u 


■ *= f A 2 / A (A ^ 0) 

(3.23) 

Sx i 

hi 

6zi 



Sx 2 

hz 

Sz2 


(3.24) 

Sx 3 

hi 

Sz 3 
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and 



*«m 

by x 

6 z x 



bu m 

6 z i 


bx { 

6y/i 


> 

H 

II® 


by 2 

bz 2 

|H 

<1 

(5x 2 

bu 2 m 

<5^2 

, A, d = 

bx 2 

by 2 

bu 2 m 


su'i 

by. 3 

bz 3 


6 x 3 

bu 3 m 

6 z 3 


bx 3 

by 3 

bu 3 m 


(3.25) 


Note that: (i) A = 0 if and only if the spatial projections of A ° , A.^, A° and Q * are 

coplanar; and (ii) For each k = 2, 3, . . . , 8, (umi)^- , (u^)q. , and (uml)Q- are defined by 
the above definition procedure except that i?, is replaced by B^ k . 

With the above preliminaries, for each m = 1,2, 3, 4, 5, ( u mx )Q. , and 

(u mz )Q- may be evaluated by 

( U mx ) Q * = («m 9 )0* = o = p (3.26) 


k=\ 


k = 1 


8 


fc=i 


Alternatively, for a flow with steep gradients or discontinuities, the simple averages in 
Eq. (3.26) may be replaced by weighted averages, i.e., 


( u mi)cj* — * 


(Umy)Q- = 


0, if0 m * = O,fc = l,2,...,8 

i /Z)t=i(^m ) ) a otherwise 

0, 

Yfk=l 0 ) ) a (Umy)Q * / Yfk=l( W ™ V Othe rwise 


(3.27a) 

if Omk = o, k - 1,2,..., 8 


(3.276) 


and 


{ u mz)Q* — 


0, if 0 m k = 0, Ar = 1, 2, . . . , 8 

J 2 k= i ( w 'm ) ) a, (‘Uml)Q- /T 2 S k=i(Wm ) ) 0 ‘ otherwise 
Here (i) a > 0 is an adjustable constant (usually a = 1 or a = 2); (ii) 


(3.27c) 


n d _ef 
v mk — 


(Umi)Q- + (Umy)Q- + (Umi)g 


.(*) 


.(*)• 


1 2 


(3.28) 


( k) • 

and (iii) for each k, Wm is the product of Q m \, 0 m2 , . . . , # m s excluding 9 m k- Note that: 
(i) to avoid dividing by zero, in practice a small positive number such as 10“ 60 is added to 
the denominators that appear in Eqs. (3.27a)-(3.27c); and (ii) Eqs. (3.27a)-(3.27c) reduce 
to Eq. (3.26) if a = 0. 
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3.5. Remarks and Discussions 

The present 3D Euler solver is formed using Eqs. (3.17) and (3.27a)-(3.27c). With 
some trivial modifications, most of the discussions about the 2D scheme given in Sec. 2.5 
and the Appendix are also applicable to the present 3D scheme. In particular, the concept 
of local and global flux conservation can also be established for the present 3D scheme by 
using a redefinition procedure similar to that presented in the Appendix. 

4. Numerical Results 


The capabilities of the present 2D and 3D schemes will be demonstrated using the 
numerical examples presented in the following subsections. 

4.1. Shock Reflection on a Flat Plate 


This steady-state test problem was proposed by Yee et al. [28]. By imposing suitable 
upstream conditions, oblique incident and reflected shocks will appear above a flat plate. 
The spatial computational domain is a 4.0 x 1.0 rectangle containing 19200 uniform rect- 
angles. For the resulting space-time mesh, (i) fi* = Q, and (ii) $2* can be divided into two 
disjoint sets Q+ and Q!_ (see Sec. 2.5). 


The flow conditions at t = 0 are [9] 

(2-9, 0.0, 1.0, 0.71428), 

' ’ v ' p ' pl \ (2.6193, -0.50632, 1.7, 1. 


ahead of the incident shock 
5282), behind the incident shock 


(4.1) 


where u, u, p and p, are x-velocity, y- velocity, mass density and static pressure, respectively. 
For t > 0, (i) the flow conditions given in the first and second rows on the right side of 
Eq. (4.1) are imposed on the left and the top boundaries, respectively; (ii) the reflecting 
boundary conditions (see the bottom half of p.124 in [9]) are imposed on the bottom 
boundary (a solid wall); and (iii) the non-reflecting conditions [9,13] are imposed on the 
right boundary (a supersonic outlet). 

Note that, for the reflecting boundary conditions used here, no mesh point lies on 
the solid wall. In addition, for each interior mesh point immediately neighboring to the 
solid wall, at the same time level there is a mirror image ghost mesh point lying just 
below the wall. Because (i) the solution values at the ghost point are assigned to be the 
mirror-image values of its corresponding interior mesh point, and (ii) one of the above two 
points belongs to 17+ while the other belongs to 171, the solution values of 17+ and 171 
are coupled by the present reflecting boundary conditions. In spite of this disadvantage, 
as explained in [9], the set of reflecting boundary conditions used here (which will also be 
used in the following numerical examples) is the most robust among several sets of the 
reflecting boundary conditions described in [9]. Note that, because the marching over 17+ 
and that over 171 are completely decoupled from each other except for the mesh points 
immediately neighboring to the solid wall, only the solution values of one of 17+ and 171 
are involved in producing Fig. 6(b), although the numerical time-marching itself involves 
both 17+ and 171. Here it should be emphasized that, for the current special problem in 
which only one straight solid wall is present, only one of 17+ and 171 needs to be used in 
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the computation if, instead, one uses the reflecting boundary conditions similar to that 
described on p.122 in [9]. 

The pressure contours generated using the present 2D scheme with a = 2 are shown in 
Fig. 6(a). The angle between the computed reflected shock and the horizontal line is 23.28°, 
which is very close to the analytical value [27]. Furthermore, as shown in Fig. 6(b), (i) the 
numerical values of the pressure coefficient at the horizontal mid-section of the rectangular 
domain agree very well with the analytical values; (ii) no numerical oscillations are detected 
near either the incident or the reflected shock; and (iii) both the incident and reflected 
shocks are resolved by a single data point. 

4.2. Shock Wave Diffraction over a Wedge 

This test problem, which was originally used by Wang [6], is based on a flow field 
given in the flow album edited by van Dyke [29]. A planar shock wave at M s = 1.3 moves 
toward a wedge with the angle 6 = 26.565 (see Fig. 7(a)). Taking advantage of symmetry, 
only half of the flow field is simulated. The spatial computational domain is a rectangle 
with —0.8 < x < 3.2 and 0 < y < 1.1, excluding the wedge. The whole domain is divided 
into 248, 750 non-uniform quadrilaterals and a = 1 is assumed. 

At t = 0. the incident planar shock is placed at x = —0.5. For t > 0, (i) the 
constant behind-the-shock flow conditions are maintained at the the left boundary; (ii) 
the reflecting boundary conditions are imposed on the upper and lower boundaries (note: 
the lower boundary is the symmetric center line), and also on the surfaces of the wedge; 
and (iii) the non-reflecting boundary conditions are imposed on the right boundary, a 
supersonic outlet. 

To enhance the visual effect, the density countours of the entire flow field at three 
different times are presented in Figs. 7(b)-(d). When the planar shock reaches the wedge, a 
circular reflection wave is generated. As the shock passes the wedge, the flow separates and 
vortices are formed around the two sharp corners. Further interaction between shocks and 
vortices produces increasingly elaborate patterns of shock waves, slip lines and vortices. 
These results agree well with the experimental result [29] except for those phenomena 
induced by the viscous effect. Here, it should be pointed out that the exact locations of 
the upper and lower walls in the experiment are not given in [29] (we only know that these 
walls are actually above and below the top and bottom edges of the photograph frame, 
respectively). As a result, the spatial domain assumed in the current simulation (which is 
slightly lareger than the photograph frame) is only an approximation of the actual physical 
domain. 

4.3. Three-Dimensional Detonation 

The 3D scheme described in Section 3 has been extended to become a solver for 
conservation laws with source terms. Previously, we have reported numerical simulations 
of ID and 2D detonation waves by using the CE/SE method [25]. Those results have been 
validated by comparing them with analytical solutions and numerical solutions reported by 
other researchers. In the present paper, 3D simulation of a detonation wave is performed 
by solving the reacting Euler equations. The chemical reactions are modeled by single-step, 
irreversible and finite-rate kinetics. Two chemical species are considered, i.e., the reactant 


and the product. The Euler equations and one species equation are solved simultaneously. 
With proper non-dimensionalization, it can he shown that the defining parameters of this 
detonation wave are the overdriven factor /, the specific heat ratio 7 , the activation energy 
E + , and the heat release rate q. I 11 the present simulation, / = 1 . 6 , 7 = 1.2, E+ = 50, 
and q = 50 are assumed. 

In the current simulation, a = 1 is assumed. Also the spatial computational domain, 
a 8 x 8 x 6 rectangular box. is divided into 6.4 million hexahedrons. Reflecting boundary 
conditions are imposed on the four lateral wall boundaries. The fresh reactant travels from 
top to bottom, and is consumed by the frame front. On the top surface, the incoming flow 
conditions are specified. On the bottom surface, a non- reflecting boundary condition is 
imposed. The coordinate system is chosen such that the frame front stays in the horizontal 
mid-section of the rectangular box. 

A snap shot of temperature countours is shown in Fig. 8 . The flow field is composed 
of the quiescent state of the reactant ahead of the shock, a flame zone with finite rate 
reaction, and the equilibrium state behind the reaction zone. Due to cellular structure of 
the detonation, the flow field is very complex. The shock front is characterized by triple 
points traveling in transverse directions. The colliding triple points create tremendous 
vortices. We observe the classical picture of “explosions within explosions” sustained by 
the propagating triple points at the detonation front. It is seen that a high-temperature 
region exists around triple points. At each collision of triple points, vortices with opposite 
signs are created and propagated downstream. Due to these vortices, unburnt reactant is 
pushed into the flame zone. The continuous burning of the pockets of the unburnt reactant 
behind the flame zone greatly extends the effective flame zone. 

5. Concluding Remarks 

In this paper, the original 2 D and 3D CE/SE Euler a-a schemes (which use triangular 
and tetrahedral meshes, respectively) were extended to solve the 2 D and 3D unsteady Euler 
equations using quadrilateral and hexahedral meshes, respectively. It has been shown that 
the present schemes retain many key advantages of other CE/SE schemes, i.e., efficient 
parallel computing, ease of implementing non-reflecting boundary conditions, high-fidelity 
solutions, and a genuinely multidimensional formulation without using Riemann solvers. 
The only key disadvantage of the present schemes (and, for that matter, any other a-a 
scheme) is that, compared with other more general CE/SE schemes such as the a-e-a-/3 
schemes [9], they allow for less freedom in adjusting numerical dissipation. As explained in 
Sec. 5.5 of [9], this inflexibility may impose a constraint on the performance of the current 
schemes in numerical simulations involving highly nonuniform meshes. 

In addition, it was pointed out that, by combining the techniques used to construct 
the present and earlier CE/SE solvers, one could easily develop 2D and 3D mixed mesh 
solvers. An advantage of using such a mixed mesh is that a geometrically complex spa- 
tial subdomain can be filled easily using triangles or tetrahedrons while a less complex 
subdomain, such as a near-wall region, can be filled using quadrilaterals or hexahedrons. 

Also, a rigorous discussion about the concept of local and global flux conservation as 
applied to the present 2 D scheme using an unstructured mesh is given in the Appendix. As 
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a part of this discussion, a post-marching procedure was introduced to handle a ‘‘solution 
decoupling" problem that may arise after a long marching involving many time steps. 
Without any exception, the discussions given in the Appendix can be extended to 3D 
easily. 


Appendix 


In this appendix, using a similar technique presented in [10], local and global flux 
conservation will be established for the present 2D scheme using an unstructured mesh. 
Also a post -marching procedure will be introduced to handle the “solution decoupling” 
problem referred to in comment (f) of Sec. 2.5. 

Note that, for the case in which the mesh decoupling referred to in comments 

(c) and (d) of Sec. 2.5 does not occur, generally the space-time computational domain 
cannot be filled by the union of a combination of nonoverlapping CCEs. As a result, 
global flux conservation cannot be established by summing over a set of local conservation 
conditions Eq. (2.13). However, even in the nondecoupling case, the computational domain 
can still be filled by the union of a combination of nonoverlapping BCEs. As a result, 
through a process of flux redefinition to be shown, one can manage to preserve the concept 
of local and global flux conservation over the BCEs and the union of any combination of 
them. 

As a preliminary, first we introduce the following definitions (see Fig. 4(b)): 

(a) For any m,£ = 1,2, 3, 4, let F^(Q*) denote the flux of h* m leaving CE(<5) through the 
top face of CE^(Q), assuming that this top face belongs to SE(Q*). Note that the top 
faces of CE^((5), (. = 1,2, 3, 4, axe the quadrilaterals A\B\QB a , A 2 B 2 QB 1 , A 3 B 3 QB 2 
and A 4 B 4 QB 3 , respectively. 

(b) For any m,l = 1,2, 3, 4, let F^A'f) denote the flux of leaving CE(Q) through 
the bottom face of CE^Q), assuming that this bottom face belongs to SE(A£*). Note 
that the bottom faces of CE^Q), t = 1,2, 3. 4, are the quadrilaterals A\B[Q' B' 4 , 
A' 2 B' 2 Q' B[, A' 2 B' 2 Q' B 2 and A' 4 B' 4 Q' B' 3 , respectively. 

(c) For any m,£ = 1,2, 3, 4 and any k = 1,2, let Fm’ e \A' e *) denote the flux of h* m 
leaving CE(Q) through its (k,£) side face, assuming that this side face belongs to 
SE(A^*). Note that the (k,£) (k = 1,2, £ — 1,2, 3, 4) side faces of CE(<3) are defined 
in Comments (e) of Sec. 2.3. 

With the above definitions, local flux conservation over CE(Q), i.e., Eq. (2.13), implies 
that 

4 

= 0 (A.l) 

e=i 


where 


s^(<3*,4*) =' fUqi + + F^(A';) + F^'\A'n 


(4.2) 


Note that Eq. (A.l) says nothing about local flux conservation over CE^(Q), — 1,2, 3, 4. 

As will be shown, local flux conservation over these BCEs can be realized with a proper 
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assignment, of "artificial fluxes' over the four interfaces that divide C'E (Q) into CEf(Q), 

? = 1 , 2 , 3 , 4 . 

To proceed, note that the boundary of each C'E t\Q) is formed by the top face, the 
bottom face and the four side faces. Among these four side faces, two are also the side 
faces of CE(Q) while the other two belong to the set of the four interfaces that divide 
CE(Q) into CEf(Q), V. = 1,2, 3, 4. Hereafter, the first pair and second pair of the above 
four side faces, respectively, are referred to as the “exterior” and “interior” side faces of 
CE e{Q). Obviously, for each m, the four terms on the right side of Eq. (A. 2) represent 
the fluxes leaving CE^Q) through its top face, bottom face, and two exterior side faces, 
respectively. 

Next, for any m = 1, 2, 3,4, let F^ 2 (Q) represent a flux (as yet to be defined explicitly) 
leaving CE L (Q) (and entering CE 2 (Q)) through the interface dividing CE^Q) and CE 2 (Q). 
F£i'\Q), F^ l (Q) and F^ l (Q) are similarly defined. In addition, for any m = 1, 2,3,4, let 


Sin(Q m , AT) + F l m 2 (Q) - F* m \Q) = 0 (A.3a) 

S 2 m (Q *, A'*) + F™(Q) - F^(Q) = 0 (A.36) 

S 3 m (Q*, AST) + Fm 4 (Q) - F 2 m \Q) = 0 (A.3c) 

S*m(Q\ A 1 ;) + F™(Q) - F 3 m 4 (Q ) = 0 (A.3d) 


Note that (i) 5^(<5*, A'*) represents the sum of the fluxes leaving CEi(<5) through its top 
face, bottom face and two exterior side faces; and (ii) Fm 2 (Q ) 311(1 -^ :1 (Q)> respectively, 
represent the fluxes leaving CEi(<5) through its two interior side faces. Thus, for each m, 
Eq. (A. 3a) represents a local flux conservation relation over CEi(Q). Similarily, for each 
m, Eqs. (A.3b)-(A.3d), represent local flux conservation relations over CE2(<5)> CE3(<3) 
and CE 4 (Q), respectively. 

Note that a summation over Eqs. (A.3a)-(A.3d) results in Eq. (A.l) — the known local 
conservation condition over CE(Q). Thus, for each m, Eqs. (A.3a)-(A.3d) contain only 
three independent conditions for four unknowns F^ 2 (Q), F^ 3 (Q), F^ 4 (Q) and F^ 1 (Q). 
In other words, there still is a degree of freedom left for these unknowns. 

To proceed, note that the interfaces that divide CE(Q) into CE<(Q), i = 1,2, 3, 4, 
all belong to SE($*). As a result, even though they are not used in the construction of 
the present scheme, the fluxes of h ^ at these interfaces can be evaluated in terms of the 
independent marching variables at point Q *. In the following discussion, the evaluated 
flux of h* m leaving CEj(Q) (and entering CE 2 (Q)) through the interface dividing CEj(Q) 
and CEo(Q) will be denoted by F^ 2 (Q*). Similarily. one also define F^ 3 (Q*), F^ 4 (Q*) 
and fT (<?•). 

With the above definitions, the degree of freedom referred to earlier is removed by 
requiring that, for each m, F^ 2 (Q), F^ 3 (Q), F^ 4 (Q) and F^ l (Q) be the solution to 
Eqs. ( A.3a)-(A.3d) with the minimal value of 

=' K 2 (<?) - oo*)] 2 + K :3 W) - F™(Q -)} 2 4 

+ IOQ) - -C(<n ] 2 + K :1 (<?> - F' m '{Q‘)Y 
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It can be shown that the last requirement amounts to imposing the extra condition 

+ OO) + + Fti\Q) 

= F'JiQ") + 0«") + Fm'iQ") + F^'(Q-) ' ' 0) 

By using Eq. (A. 5) and any three of Eqs. (A.3a)-(A.3d), for each m, F^ 2 (Q), F^ 3 (Q), 
F'm 4 (Q) and F^ l (Q) can be uniquely defined in terms of known parameters F^ 2 {Q*), 
F™(Q*), F 3 m \Q*), F^(Q’), and S^( Q’ . A'f). £ = 1,2, 3, 4. 

Next, note that a space-time region may be the common BCE of two different mesh 
points (these two mesh points are referred to as the cohosts of the common BCE). As an 
example, the space-time cylinder A i B\QB i A' l B[ Q' B^ depicted in Fig. 4(b) was designated 
as CE\(Q). However, it also can be designated as a BCE of point A\, say CEi(Ai). As will 
be shown in the following remarks, for each m, how the flux is assigned to each face of the 
space-time cylinder, along with the resulting flux conservation relation over the cylinder, 
is dependent on whether it is designated as CEx(Q) or CEi(Ai): 

(a) At the top face of CEi(Q) (CEi(Ai)), the flux is evaluated assuming that the face 
belongs to SE (Q*) (SE (4*)). 

(b) At the bottom face of CEi(<5) (CEj(Ai)), the flux is evaluated assuming that the face 
belongs to SE(A;*) (SE(Q'*)) 

(c) The exterior (interior) side faces of CEi(Q) are the interior (exterior) side faces of 

CE 1 (A 1 ). 

(d) At each of the exterior side faces of CE^Q) (CEi(Ai)), the flux is evaluated assuming 
that the side face belongs to SE(Ai*) (SE(Q*)). 

(e) A local conservation condition over CEi(Aj) (different from that over CEj(Q), i.e., 
Eq. (A. 3a)) will result if the artificial flux at each interior side face of CE^Aj) is also 
assigned using a procedure parallel to that used to assign the flux at each interior side 
face of CEj((5). 

Consider a common BCE of two cohosts lying in the interior of the computational 
domain. From the above discussion, one concludes that, for each m, (i) two different fluxes 
are assigned to each face of the BCE, and (ii) corresponding to the two cohosts, there are 
two different conservation relations over this BCE. Hereafter, the simple average of the 
two fluxes at each face will be referred to as the generalized flux at this face. By summing 
the two local conservation relations over the BCE, one concludes that the total generalized 
flux leaving the BCE through its boundary vanishes. 

Furthermore, note that: (i) only one generalized flux is defined at any interface divid- 
ing two neighboring BCEs; and (ii) the generalized flux leaving a BCE through an interface 
dividing this BCE and a neighboring BCE is the negative of the generalized flux leaving 
the neighboring BCE through the same interface. Thus, one arrives at the following global 
flux conservation relation: for each m, the total generalized flux leaving the boundary of 
any space-time region that is the union of any combination of BCEs (with each of these 
BCEs having two interior cohosts) vanishes. 

To proceed further, note that, for each m, corresponding to its two cohosts, the 
boundary of a BCE is assigned two sets of fluxes. Because of the space-time staggering 
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nature of the stencil of the present scheme, the above two sets along with the solution values 
at its two cohosts may become decoupled locally after many marching steps. Note that one 
may argue that this decoupling does not matter, because the amount, of decoupling usually 
is of the order of the discrepancy between the numerical solution and the exact solution 
and, as such, it does not exacerbate the actual simulation errors. However, in practice, 
the decoupling can cause a substantial problem in solution display. The decoupling can 
manifest itself as what appear to be small- wavelength oscillations when the solution at 
the final time level is displayed using the solution values of both and D*_. As will 
be shown immediately, not only does the above definition of a unique generalized flux at 
any boundary of a BCE provide a way to avoid the problem of “flux decoupling”, it also 
provides a way to handle the problem of “solution decoupling”. 

Consider the top face of any BCE with two cohosts. For any m, the two fluxes assigned 
to this face, respectively, are evaluated assuming that the face belongs to the SEs of its 
two cohosts, respectively. It can be shown that these two fluxes, respectively, are equal 
to the area of the face multiplied by the two values of u m at the centroid of the top face 
evaluated assuming that the centroid belongs to the two cohosts, respectively. Let the 
simple average of the above two values of u m be referred to as the coupled solution value 
of u m at the centroid of the top face of this BCE. Then it can easily be shown that, for 
each m, the generalized flux at this face is simply the area of the face multiplied by the 
new solution value. Also, because of how they are defined, solution decoupling generally 
is no longer a problem if the numerical data are taken from these new solution values. 

Finally, it should be emphasized that the above definition of generalized fluxes and 
coupled solution values, by no means implies any change in the marching scheme. In fact, 
evaluation of the locations of the centroids of the top faces of the BCEs along with that 
of the associated coupled solution values represents only a post-marching procedure. 
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Figure 1 . — A surface element on the boundary S (V) 
of a volume V in a space-time E 2 . 
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Figure 2. — Space-time geometry of the conventional 
finite volume method in £ 2 - (a) A rectangle in E 2 . 

(b) A spatial cylinder aligned in the x-direction, 

(c) A regular space-time mesh. 
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Figure 3. — The SEs and CEs of the a scheme, (a) A staggered 
space-time mesh, (b) SE(j,n). (c) CE_(j,n). (d) CE + (j,n). (e) CE(j,n). 
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Figure 4. — Space-time geometry of the 2D scheme. 

(a) Representative grid points in the x-y plane. 

(b) SEs and CEs. (c) Spatial translation of the 
quadrilateral A] A£ A 3 A 4 . 
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Figure 5. — Representative grid points 
in the x-y-z space. 
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Fig. 6: The Euler solution of a steady-state shock reflection problem: (a) pressure contours; (b) 
pressure coefficient distribution at the mid-section of the computation domain (y=0.5). 
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Fig. 7: Schematic and density contours at three different times compared with the experimental 
photographs. 







